"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBTools.py 990 2012-02-15 13:52:11Z bgoli $)

"""

import os, time, re, cPickle, copy
cDir = os.path.dirname(os.path.abspath(os.sys.argv[0]))
import numpy
import CBModel
from CBDataStruct import StructMatrixLP # legacy
from CBConfig import __CBCONFIG__ as __CBCONFIG__
__DEBUG__ = __CBCONFIG__['DEBUG']
__version__ = __CBCONFIG__['VERSION']

def storeObj(obj, filename):
    """
    Stores a Python *obj* as a serialised binary object in *filename.dat*

    """
    filename = filename+'.dat'
    F = file(filename, 'wb')
    cPickle.dump(obj, F, protocol=2)
    print 'Object serialised as %s' % filename
    F.close()

def loadObj(filename):
    """
    Loads a serialised Python cPickle from *filename.dat* returns the Python object(s)
    """
    return loadObject(filename)

def loadObject(filename):
    """
    Loads a serialised Python cPickle from *filename.dat* returns the Python object(s)
    
    **Note** this is being deprecated and refactored ... use loadObj()
    """    
    assert os.path.exists(filename), '\nTry again mate!'
    F = file(filename, 'rb')
    obj = cPickle.load(F)
    F.close()
    return obj

def addStoichToFBAModel(fm):
    """
    Build stoichiometry: this method has been refactored into the model class - cmod.buildStoichMatrix()
    """
    fm.buildStoichMatrix()



def addSinkReaction(fbam, species, lb=0.0, ub=1000.0):
    """
    Adds a sink reactions that consumes a model *species* so that X -->

     - *fbam* an fba model object
     - *species* a valid species name
     - *lb* lower flux bound [default = 0.0]
     - *ub* upper flux bound [default = 1000.0]

    """
    assert species in fbam.getSpeciesIds(), '\n%s is not a valid species' % species
    if lb < 0.0:
        reversible = True
    else:
        reversible = False
    Rname = species+'_sink'
    R = CBModel.Reaction(Rname, name='%s sink reaction' % species, reversible=reversible)
    Su = CBModel.Reagent(Rname+species, species, -1.0)
    R.addReagent(Su)
    R.is_exchange = True

    clb = CBModel.FluxBound(Rname+'_lb', Rname, 'greaterEqual', lb)
    cub = CBModel.FluxBound(Rname+'_ub', Rname, 'lessEqual', ub)

    fbam.addReaction(R)
    fbam.addFluxBound(clb)
    fbam.addFluxBound(cub)

    print '\n***\nCreated new reaction %s with bounds (%s : %s)\n***\n' % (Rname, lb, ub)


# TODO: check this
def addSourceReaction(fbam, species, lb=0.0, ub=1000.0):
    """
    Adds a source reactions that produces a model *species* so that --> X

     - *fbam* an fba model object
     - *species* a valid species name
     - *lb* lower flux bound [default = 0.0]
     - *ub* upper flux bound [default = 1000.0]
     
    Note reversiblity is determined by the lower bound, default 0 = irreversible. If
    negative then reversible.    

    """
    assert species in fbam.getSpeciesIds(), '\n%s is not a valid species' % species
    if lb < 0.0:
        reversible = True
    else:
        reversible = False
    Rname = species+'_src'
    R = CBModel.Reaction(Rname, name='%s source reaction' % species, reversible=reversible)
    Su = CBModel.Reagent(Rname+species, species, 1.0)
    R.addReagent(Su)
    R.is_exchange = True

    clb = CBModel.FluxBound(Rname+'_lb', Rname, 'greaterEqual', lb)
    cub = CBModel.FluxBound(Rname+'_ub', Rname, 'lessEqual', ub)

    fbam.addReaction(R)
    fbam.addFluxBound(clb)
    fbam.addFluxBound(cub)

    print '\n***\nCreated new reaction %s with bounds (%s : %s)\n***\n' % (Rname, lb, ub)



def removeOphanedReactions(fbam):
    orphaned_list = []
    for rr in range(fbam.N.array.shape[0]):
        if (fbam.N.array[rr, :] != 0.0).sum() == 1:
            if __DEBUG__: print fbam.N.array[rr,:]
            if __DEBUG__: print fbam.N.row[rr]
            for c in range(len(fbam.N.array[rr,:])):
                if fbam.N.array[rr,c] != 0.0:
                    orphaned_list.append(fbam.N.col[c])
                    if __DEBUG__: print '\t', fbam.N.col[c]
            fbam.N.array[rr,:] = 0.0
    return orphaned_list

def FixReversibility(fbam, auto_correct=False):
    """
    Set fluxbound lower bound from reactions reversibility information.

     - *fbam* and FBAModel instance
     - *auto_correct* (default=False) if True automatically sets lower bound to zero if required, otherwise prints a warning if false.

    """
    for c in fbam.flux_bounds:
        R = c.reaction
        ##  print R
        O = c.operation
        ##  print O
        V = c.value
        ##  print V
        R_obj = fbam.reactions[fbam.getReactionIds().index(c.reaction)]
        RE = R_obj.reversible
        ##  print RE
        if O in ['greater', 'greaterEqual']:
            if not RE and float(V) < 0.0:
                print 'Warning %s is not reversible and lower bound is %s.' % (R, V)
                if auto_correct:
                    print 'Resetting %s lower bound (%s) to zero' % (R, V)
                    c.value = '0.0'
                else:
                    print 'Reaction (%s) reversible=%s inconsistent with fluxbound lower bound (%s) run with auto_correct=True to reset lower bound.' % (R, RE, V)
                    time.sleep(1)


def SplitReversibleReactions():
    """
    Split a (set of) reactions into reversible reactions.

    R1: A = B
    R1f: A -> B
    R1r: B -> A

    """
    raise NotImplementedError

def getReversible():
    """
    List reversible etc.

    """
    raise NotImplementedError

def getReducedCosts(fba):
    """
    Get a dictionary of reduced costs for each reaction/flux

    """
    output = {}
    for r in fba.reactions:
        output.update({r.getPid() : r.reduced_cost})
    return output

def setReducedCosts(fba, reduced_costs):
    """
    For each reaction/flux, sets the attribute "reduced_cost" from a dictionary of
    reduced costs

     - *fba* an fba object
     - *reduced_costs* a dictionary of {reaction : value} pairs

    """
    for r in fba.reactions:
        if reduced_costs.has_key(r.getPid()):
            r.reduced_cost = reduced_costs[r.getPid()]
        else:
            r.reduced_cost = None

def exportLabelledArray(arr, names, fname, sep=',', format='%f'):
    if names != None:
        assert arr.shape[0] == len(names), '\n ...  rows must equal number of names!'
    F = file(fname, 'w')
    cntr = 0
    for r in range(arr.shape[0]):
        if names != None:
            F.write(('%s'+sep) % names[r])
        for c in range(arr.shape[1]):
            if c < arr.shape[1]-1:
                F.write((format+sep) % arr[r,c])
            else:
                F.write((format+'\n') % arr[r,c])
        cntr += 1
        if cntr >= 250:
            F.flush()
            cntr = 1
    F.write('\n')
    F.flush()
    F.close()
    print 'exported to %s' % fname

def exportLabelledArrayWithHeader(arr, names=None, header=None, fname=None, sep=',', format='%f'):
    """
    Export an array with names=['row','names'] and header=['column','names']
    """
    if fname == None:
        raise RuntimeError, '\nPlease enter a filename'
    if names != None:
        assert arr.shape[0] == len(names), '\n ...  rows must equal number of names!'
        
    if header != None:
        assert arr.shape[1] == len(header), '\n ...  cols must equal number of header names!'
    F = file(fname, 'w')
    cntr = 0
    if header != None:
        if names != None:
            hstr = ' '+sep
        else:
            hstr = ''
        for h in header:
            hstr += str(h)+sep
        hstr = hstr[:-1]+'\n'
        F.write(hstr)
        del hstr
    for r in range(arr.shape[0]):
        if names != None:
            F.write(('%s'+sep) % names[r])
        for c in range(arr.shape[1]):
            if c < arr.shape[1]-1:
                F.write((format+sep) % arr[r,c])
            else:
                F.write((format+'\n') % arr[r,c])
        cntr += 1
        if cntr >= 250:
            F.flush()
            cntr = 1
    F.write('\n')
    F.flush()
    F.close()
    print 'exported to %s' % fname

def exportLabelledArray_altzero(arr, names, fname, sep=',', format='%f'):
    if names != None:
        assert arr.shape[0] == len(names), '\n ...  rows must equal number of names!'
    F = file(fname, 'w')
    cntr = 0
    for r in range(arr.shape[0]):
        if names != None:
            F.write(('%s'+sep) % names[r])
        for c in range(arr.shape[1]):
            if c < arr.shape[1]-1:
                if arr[r,c] == 0.0:
                    F.write('0.0'+sep)
                else:
                    F.write((format+sep) % arr[r,c])
            else:
                if arr[r,c] == 0.0:
                    F.write('0.0\n')
                else:
                    F.write((format+'\n') % arr[r,c])
        cntr += 1
        if cntr >= 250:
            F.flush()
            cntr = 1
    F.write('\n')
    F.flush()
    F.close()
    print 'exported to %s' % fname

def exportLabelledLinkedList(arr, names, fname, sep=',', format='%s', appendlist=False):
    """
    Write a 2D linked list [[...],[...],[...],[...]] and optionally a list of row labels to file:
    
     - *arr* the linked list
     - *names* the list of row names
     - *fname* the output filename
     - *sep* [default=','] the column separator
     - *format* [default='%s'] the output number format
     - *appendlist* [default=False] if True append the array to *fname* otherwise create a new file
     
    """
    if names != None:
        assert len(arr) == len(names), '\n ...  rows must equal number of names!'
    if not appendlist:
        F = file(fname, 'w')
    else:
        F = file(fname, 'a')
    cntr = 0
    for r in range(len(arr)):
        if names != None:
            F.write(('%s'+sep) % names[r])
        col_l = len(arr[0])
        for c in range(col_l):
            if c < col_l-1:
                if arr[r][c] == 0.0:
                    F.write('0.0'+sep)
                else:
                    F.write((format+sep) % arr[r][c])
            else:
                if arr[r][c] == 0.0:
                    F.write('0.0\n')
                else:
                    F.write((format+'\n') % arr[r][c])
        cntr += 1
        if cntr >= 250:
            F.flush()
            cntr = 1
    ##  F.write('\n')
    F.flush()
    F.close()
    del arr
    if not appendlist:
        print 'exported to %s' % fname

def exportLabelledArrayWithHeader2CSV(arr, names=None, header=None, fname=None):
    if fname == None:
        raise RuntimeError, '\nPlease enter a filename'    
    fname += '.csv'
    exportLabelledArrayWithHeader(arr, names, header, fname, sep=',', format='%f')

def exportLabelledArray2CSV(arr, names, fname):
    fname += '.csv'
    exportLabelledArray(arr, names, fname, sep=',', format='%f')

def exportArray2CSV(arr, fname):
    fname += '.csv'
    exportLabelledArray(arr, None, fname, sep=',', format='%f')

class ComboGen:
    combo = None
    combo_int = None

    def __init__(self):
        self.combo = []

    def addCombination(self,data):
        self.combo.append(','.join(data))

    def uniqueCombinations(self,data,number,temp=[]):
        if not number:
            self.addCombination(temp)
        for i in range(len(data)):
            temp.append(data[i])
            self.uniqueCombinations(data[i+1:],number-1,temp)
            temp.pop()

    def numberifyComb2Int(self):
        tmp = []
        for c in self.combo:
            tmp.append(tuple([int(c) for c in c.split(',')]))
        self.combo_int = tmp

def stringReplace(fbamod, old, new, target):
    """
    This is alpha stuff, target can be:

     - 'species'
     - 'reactions'
     - 'constraints'
     - 'objectives'
     - 'all'

    """
    print 'stringReplace is relatively new and UNTESTED'
    fbamod.id = fbamod.id.replace(old, new)
    if target == 'species' or target == 'all':
        for s in fbamod.species:
            s.id = s.id.replace(old, new)
    if target == 'reactions' or target == 'all':
        for s in fbamod.reactions:
            s.id = s.id.replace(old, new)
            for r in s.reagents:
                r.id = r.id.replace(old, new)
    if target == 'constraints' or target == 'all':
        for s in fbamod.flux_bounds:
            s.id = s.id.replace(old, new)
            s.reaction = s.reaction.replace(old, new)
    if target == 'objectives' or target == 'all':
        for s in fbamod.objectives:
            s.id = s.id.replace(old, new)
            for f in s.fluxObjectives:
                f.id = f.id.replace(old, new)
                f.reaction = f.reaction.replace(old, new)
    return fbamod

def getBoundsDict(fbamod, hasstring=None):
    """
    Return a dictionary of reactions&bounds
    """
    rBdic = {}
    for r in fbamod.getReactionIds(hasstring=hasstring):
        name, lb, ub, eq = fbamod.getReactionBounds(r)
        rBdic.update({name : {'lb' : lb,
                              'ub' : ub,
                              'eq' : eq
                             }
                    })
    return rBdic

def getExchBoundsDict(fbamod):
    """
    Return a dictionary of all exchange reactions (as determined by the is_exchange attribute of Reaction)
    """
    rBdic = {}
    for r in fbamod.getReactionIds(hasstring=None):
        name, lb, ub, eq = fbamod.getReactionBounds(r)
        rBdic.update({name : {'lb' : lb,
                              'ub' : ub,
                              'eq' : eq
                             }
                    })
    for r in fbamod.reactions:
        if not r.is_exchange:
            rBdic.pop(r.getPid())
    return rBdic

def processBiGGchemFormula(fba):
    """
    Disambiguates the overloaded BiGG name NAME_CHEMFORMULA into

     - *species.name* NAME
     - *species.chemFormula* CHEMFORMULA

    """
    for s in fba.species:
        try:
            tmp = s.name
            tmp2 = tmp.split('_')
            tmp2.reverse()
            CF = tmp2[0]
            NM = tmp.replace('_%s' % CF, '')
            if __DEBUG__:
                print NM, CF
        except Exception, ex:
            print ex
            time.sleep(1)
            tmp = None
            NM = s.name
            CF = ''
        del tmp, tmp2
        s.chemFormula = CF.strip()
        s.name = NM.strip()

def processBiGGannotationNote(fba, annotation_key='note'):
    """
    Parse the HTML formatted reaction information stored in the BiGG notes field.
    This function is being deprecated and replaced by `CBTools.processSBMLAnnotationNotes()`

     - requires an *annotation_key* which contains a BiGG HTML fragment

    """
    print '\nDeprecation warning:\nCBTools.processBiGGannotationNote() is being replaced with CBTools.processSBMLAnnotationNotes'
    
    html_p = re.compile("<html:p>.*?</html:p>")
    for r in fba.reactions:
        new_ann = {}
        if r.annotation.has_key(annotation_key):
            hPs = re.findall(html_p, r.annotation.pop(annotation_key))
            if __DEBUG__: print hPs
            for p in hPs:
                ps = p.replace('<html:p>','').replace('</html:p>','').replace('&lt;','<').replace('&gt;','>').split(':',1)
                if len(ps) == 2:
                    new_ann.update({ps[0].strip() : ps[1].strip()})
            r.annotation.update(new_ann)
        if __DEBUG__: print r.annotation

def processSBMLAnnotationNotes(fba, annotation_key='note'):
    """
    Parse the HTML formatted reaction information stored in the SBML notes field currently
    processes BiGG and PySCeSCBM style annotations it looks for the the annotation indexed
    with the *annotation_key* 

     - *annotation_key* [default='note'] which contains a HTML/XHTML fragment in BiGG/PySCeSCBM format

    """
    html_p = re.compile("<p>.*?</p>")
    html_span = re.compile("<span>.*?</span>")
    html_bigg_p = re.compile("<html:p>.*?</html:p>")
    for r in fba.reactions:
        if r.annotation.has_key(annotation_key):
            new_ann = {}
            if '<span xmlns="http://www.w3.org/1999/xhtml">' in r.annotation[annotation_key]:
                hPs = re.findall(html_p, r.annotation.pop(annotation_key).replace('\n',''))
                if __DEBUG__: print hPs
                for p in hPs:
                    ps = re.findall(html_span, p)
                    ps = [p.replace('<span>','').replace('</span>','').replace('&lt;','<').replace('&gt;','>').strip() for p in ps]
                    if len(ps) == 2:
                        new_ann.update({ps[0] : ps[1]})
            else:
                hPs = re.findall(html_bigg_p, r.annotation.pop(annotation_key))
                if __DEBUG__: print hPs
                for p in hPs:
                    ps = p.replace('<html:p>','').replace('</html:p>','').replace('&lt;','<').replace('&gt;','>').split(':',1)
                    if len(ps) == 2:
                        new_ann.update({ps[0].strip() : ps[1].strip()})                
            r.annotation.update(new_ann)
            if __DEBUG__: print r.annotation
    for s in fba.species:
        if s.annotation.has_key(annotation_key):
            new_ann = {}
            if '<span xmlns="http://www.w3.org/1999/xhtml">' in s.annotation[annotation_key]:
                hPs = re.findall(html_p, s.annotation.pop(annotation_key).replace('\n',''))
                if __DEBUG__: print hPs
                for p in hPs:
                    ps = re.findall(html_span, p)
                    ps = [p.replace('<span>','').replace('</span>','').replace('&lt;','<').replace('&gt;','>').strip() for p in ps]
                    if len(ps) == 2:
                        new_ann.update({ps[0] : ps[1]})
            else:
                hPs = re.findall(html_bigg_p, s.annotation.pop(annotation_key))
                if __DEBUG__: print hPs
                for p in hPs:
                    ps = p.replace('<html:p>','').replace('</html:p>','').replace('&lt;','<').replace('&gt;','>').split(':',1)
                    if len(ps) == 2:
                        new_ann.update({ps[0].strip() : ps[1].strip()})    
            s.annotation.update(new_ann)
        if s.annotation.has_key('chemFormula') and (s.chemFormula == None or s.chemFormula == ''):
            s.chemFormula = s.annotation.pop('chemFormula')
            if __DEBUG__: print r.annotation        

def processExchangeReactions(fba, key):
    """
    Extract exchange reactions from model using *key* and return:

     - a dictionary of all exchange reactions without *medium* reactions
     - a dictionary of *medium* exchange reactions (negative lower bound)

    """
    # extract all exchange bounds
    if key == None:
        fexDic = getExchBoundsDict(fba)
    else:
        fexDic = getBoundsDict(fba, hasstring=key)


    # extract the medium (exchange fluxes that allow uptake)
    MediumAll = []
    Medium = []
    for r in fexDic:
        if fexDic[r]['lb'] < 0.0:
            MediumAll.append((r, fexDic[r]['lb'], fexDic[r]['ub']))
            Medium.append(r)
        if __DEBUG__: print r, fexDic[r]
    # remove medium from bounds dictionary and place in medium dict
    mediumDic = {}
    for m in Medium:
        mediumDic.update({m : fexDic.pop(m)})
    if __DEBUG__:
        print '\nMedium'
        for m in MediumAll:
            print m
        print 'mediumDic'
        print mediumDic
        print '\nr in fexDic'
        for r in mediumDic.keys():
            print r, fexDic.has_key(r)
    return fexDic, mediumDic

def runInputScan(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False):
    raise RuntimeError, 'This method has been refactored to solver specific modules, access via CBSolve'

def runInputScanWithFVA(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False, NaNonFail=False):
    raise RuntimeError, 'This method has been refactored to solver specific modules, access via CBSolve'


def generateInputScanReports(fba, exDict, mediumDict, optimal_growth_rates, wDir, tag=''):
    modName = fba.inputfile_id
    modName += tag
    rnames = fba.getReactionNames()
    rid = fba.getReactionIds()
    F = file(os.path.join(wDir, '%s.medium.csv' % modName), 'w')
    F.write('J, lb, ub\n')
    for r in mediumDict.keys():
        RN = rnames[rid.index(r)]
        F.write('%s, %s, %s, %s\n' % (r, mediumDict[r]['lb'], mediumDict[r]['ub'], RN))
    F.write('\n')
    F.close()

    F = file(os.path.join(wDir, '%s.exchange.csv' % modName), 'w')
    F.write('J, lb, ub\n')
    for r in exDict.keys():
        RN = rnames[rid.index(r)]
        F.write('%s, %s, %s, %s\n' % (r, exDict[r]['lb'], exDict[r]['ub'], RN))
    F.write('\n')
    F.close()

    F = file(os.path.join(wDir, '%s.optima.csv' % modName), 'w')
    F.write('J, lb, ub, optimum, "name"\n')
    if __DEBUG__:
        print rnames
        print rid
    for r in optimal_growth_rates.keys():
        RN = rnames[rid.index(r)]
        F.write('%s, %s, %s, %s, "%s"\n' % (r, exDict[r]['lb'], exDict[r]['ub'], optimal_growth_rates[r], RN))
    F.write('\n')
    F.close()

def getAllReactionsAssociatedWithGene(fba, gene, gene_annotation_key='GENE ASSOCIATION'):
    out = []
    for r in fba.reactions:
        GA = None
        if r.annotation.has_key(gene_annotation_key):
            GA = gene_annotation_key
        elif r.annotation.has_key('GENE ASSOCIATION'):
            GA = 'GENE ASSOCIATION'
        elif r.annotation.has_key('GENE_ASSOCIATION'):
            GA = 'GENE_ASSOCIATION'
        if GA != None:
            if gene in r.annotation[GA]:
                out.append(r.getPid())
    return out

def scanForReactionDuplicates(f, ignore_coefficients=False):
    """
    This method uses uses a brute force apprach to finding reactions with matching
    stoichiometry
    """
    duplicates = []
    for r in f.reactions:
        Rref = r.getReagentRefs()
        Rref.sort()
        refspecies = ''
        for s in Rref:
            refspecies += '%s:' % s
        refspecies = refspecies[:-1]
        for r2 in f.reactions:
            Rtest = r2.getReagentRefs()
            Rtest.sort()
            if Rref == Rtest and r.id != r2.id:
                if not ignore_coefficients:
                    go = True
                    for rgid in Rtest:
                        if float(r.getReagentWithSpeciesRef(rgid).coefficient) != float(r2.getReagentWithSpeciesRef(rgid).coefficient):
                            go = False
                            break
                    if go:
                        dup = [r.id, r2.id,]
                        dup.sort()
                        dup = dup+[refspecies, f.getReaction(dup[0]).getName(), f.getReaction(dup[1]).getName()]
                        if dup not in duplicates:
                            duplicates.append(dup)
                else:
                    dup = [r.id, r2.id,]
                    dup.sort()
                    dup = dup+[refspecies, f.getReaction(dup[0]).getName(), f.getReaction(dup[1]).getName()]
                    if dup not in duplicates:
                        duplicates.append(dup)
    for d in duplicates:
        print(d)
    print('\nFound %s pairs of duplicate reactions' % len(duplicates))        
    return duplicates

def countedPause(Tsec):
    print '\nPausing ...   ',
    for t in range(Tsec,-1,-1):
        print '\b\b\b%2i' % t,
        time.sleep(1)
    print '\b\b%s' % 'done.'


# these methods will eventually go into CBDbase and/or CBNettools modules


# these methods will eventually go into CBDbase and/or CBNettools modules



def addProteinCostAnnotation(F, reaction_gene_map, gene_peptide_map, gene_annotation_key='GENE ASSOCIATION'):
    reaction_gene_map = reaction_gene_map.copy()
    gene_peptide_map = gene_peptide_map.copy()
    for r in F.reactions:
        if reaction_gene_map[r.getPid()] != None:
            r.annotation.update({'CBM_PEPTIDE_LENGTHS' : [(g, gene_peptide_map[g]) for g in reaction_gene_map[r.getPid()]]})
        else:
            r.annotation.update({'CBM_PEPTIDE_LENGTHS' : []})
        GA = None
        if r.annotation.has_key(gene_annotation_key):
            GA = gene_annotation_key
        elif r.annotation.has_key('GENE ASSOCIATION'):
            GA = 'GENE ASSOCIATION'
        elif r.annotation.has_key('GENE_ASSOCIATION'):
            GA = 'GENE_ASSOCIATION'
        if GA != None:
            newAss = r.annotation[GA]
            for rl in r.annotation['CBM_PEPTIDE_LENGTHS']:
                ##  print rl
                ##  print newAss
                newAss = newAss.replace(rl[0], str(rl[1]))
            newAss = newAss.replace('and','+').replace('(','').replace(')','')
            r.annotation.update({'CBM_PEPTIDE_SUBSTITUTED_LENGTHS' : newAss})
            ##  print r.annotation
            if 'or' in r.annotation[GA]:
                ##  print ' '
                ##  print r.annotation[GA]
                ##  print r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS']
                ##  print r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS'].split('or')
                splitor = []
                for ass in r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS'].split('or'):
                    exec('splitor.append(%s)' % ass)
                ##  print splitor
                r.annotation.update({'CBM_PEPTIDE_LENGTH_MIN' : min(splitor)})
                r.annotation.update({'CBM_PEPTIDE_LENGTH_MAX' : max(splitor)})
            else:
                if r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS'] != None and r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS'] != '':
                    maxL = 0
                    ##  print r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS']
                    exec('maxL = %s' % r.annotation['CBM_PEPTIDE_SUBSTITUTED_LENGTHS'])
                    r.annotation.update({'CBM_PEPTIDE_LENGTH_MIN' : maxL})
                    r.annotation.update({'CBM_PEPTIDE_LENGTH_MAX' : maxL})
                else:
                    r.annotation.update({'CBM_PEPTIDE_LENGTH_MIN' : None})
                    r.annotation.update({'CBM_PEPTIDE_LENGTH_MAX' : None})
            ##  print r.getPid(), r.annotation['CBM_PEPTIDE_LENGTH_MIN'], r.annotation['CBM_PEPTIDE_LENGTH_MAX']
        else:
            r.annotation.update({'CBM_PEPTIDE_SUBSTITUTED_LENGTHS' : None})
            r.annotation.update({'CBM_PEPTIDE_LENGTH_MIN' : None})
            r.annotation.update({'CBM_PEPTIDE_LENGTH_MAX' : None})

def scaleProteinCostAnnotation(fba):
    LEN_MAX = max([r.annotation['CBM_PEPTIDE_LENGTH_MAX'] for r in fba.reactions if r.annotation['CBM_PEPTIDE_LENGTH_MAX'] != None])
    LEN_AVG = numpy.average([r.annotation['CBM_PEPTIDE_LENGTH_MAX'] for r in fba.reactions if r.annotation['CBM_PEPTIDE_LENGTH_MAX'] != None])

    # cost function (AAlen/avg(AAlen) or 1.0
    for R in fba.reactions:
        if R.annotation['CBM_PEPTIDE_LENGTH_MAX'] == None:
            R.annotation['CBM_PEPTIDE_COST'] = 1.0
        else:
            R.annotation['CBM_PEPTIDE_COST'] = R.annotation['CBM_PEPTIDE_LENGTH_MAX']/LEN_AVG
        R.annotation['CBM_AVG_PEPTIDE_LENGTH'] = LEN_AVG
    print '\nAverage protein length:', LEN_AVG
        ##  print R.getPid(), R.annotation['CBM_PEPTIDE_COST'], R.annotation['CBM_PEPTIDE_LENGTH_MAX']

# Standard amino acid abbreviations with ATP formation/production cost: Stouthamer (1973) Table 4 page 555
amino_acid_table = [['A', 'Ala', 'Alanine', '+1'],
                    ['R', 'Arg', 'Arginine', '-3'],
                    ['N', 'Asn', 'Asparagine', '-2'],
                    ['D', 'Asp', 'Aspartic acid', '+0'],
                    ['C', 'Cys', 'Cysteine', '-3'],
                    ['E', 'Glu', 'Glutamic acid', '+1'],
                    ['Q', 'Gln', 'Glutamine', '+0'],
                    ['G', 'Gly', 'Glycine', '+0'],
                    ['H', 'His', 'Histidine', '-7'],
                    ['I', 'Ile', 'Isoleucine', '-1'],
                    ['L', 'Leu', 'Leucine', '+3'],
                    ['K', 'Lys', 'Lysine', '+0'],
                    ['M', 'Met', 'Methionine', '-4'],
                    ['F', 'Phe', 'Phenylalanine', '-2'],
                    ['P', 'Pro', 'Proline', '+0'],
                    ['S', 'Ser', 'Serine', '+0'],
                    ['T', 'Thr', 'Threonine', '-2'],
                    ['W', 'Trp', 'Tryptophan', '-5'],
                    ['Y', 'Tyr', 'Tyrosine', '-2'],
                    ['V', 'Val', 'Valine', '+2']
                   ]




def getModelGenesPerReaction(fba, gene_pattern='(b\w*\W)', gene_annotation_key='GENE ASSOCIATION'):
    react_gene = {}
    gene_re = re.compile(gene_pattern)

    for r in fba.reactions:
        GA = None
        if r.annotation.has_key(gene_annotation_key):
            GA = gene_annotation_key
        elif r.annotation.has_key('GENE ASSOCIATION'):
            GA = 'GENE ASSOCIATION'
        elif r.annotation.has_key('GENE_ASSOCIATION'):
            GA = 'GENE_ASSOCIATION'
        genes = re.findall(gene_re, r.annotation[GA])
        genes = [g.replace('(','').replace(')','').strip() for g in genes]
        ##  print r.annotation['GENE ASSOCIATION']
        if len(genes) == 0:
            ##  print '\n'
            ##  print 'GA:', r.annotation['GENE ASSOCIATION']
            ##  print r.getPid(), genes
            ##  raw_input('x')
            genes = None
        ##  print r.getPid(), genes
        ##  raw_input()
        react_gene.update({r.getPid() : genes})
    return react_gene

def getReactionsPerGene(react_gene):
    gene_react = {}
    no_gene = []
    for R in react_gene:
        if react_gene[R] == None:
            print 'Reaction %s has no gene associated with it' % R
            no_gene.append(R)
        else:
            for G in react_gene[R]:
                ##  GK = G.replace('(','').replace(')','').strip()
                if gene_react.has_key(G):
                    print 'Updating gene %s with reaction %s' % (G, R)
                    gene_react[G].append(R)
                else:
                    print 'Adding gene %s to gene_react' % G
                    gene_react.update({G : [R]})
    genes = gene_react.keys()
    return gene_react, genes, no_gene

def removeFixedSpeciesReactions(f):
    """
    This function is a hack that removes reactions which only have boundary species as reactants
    and products. These are typically gene associations encoded in the Manchester style and there
    is probably a better way of working around this problem ...

     - *f* an instantiated fba model object

    """
    c_react = []
    for re in f.reactions:
        lsa = numpy.array([f.getSpecies(r.species_ref).is_boundary for r in re.reagents])
        if lsa.all():
            c_react.append(re.getPid())
    for r in c_react:
        f.deleteReactionAndBounds(r)


def addFluxAsActiveObjective(f, reaction_id, osense, coefficient=1):
    """
    Adds a flux as an active objective function

     - *reaction_id* a string containing a reaction id
     - *osense* objective sense must be **maximize** or **minimize**
     - *coefficient* the objective funtion coefficient [default=1]

    """
    osense = osense.lower()
    if osense == 'max': osense = 'maximize'
    if osense == 'min': osense = 'minimize'    
    if osense in ['maximise', 'minimise']:
        osense = osense.replace('se','ze')
    assert osense in ['maximize', 'minimize'], "\nosense must be ['maximize', 'minimize'] not %s" % osense   
    assert reaction_id in [r.getPid() for r in f.reactions], '\n%s is not avalid reaction' % reaction_id
    n_obj = CBModel.Objective(reaction_id+'_objf',osense)
    n_flux_obj = CBModel.FluxObjective(reaction_id+'_fluxobj', reaction_id, coefficient)
    n_obj.addFluxObjective(n_flux_obj)
    f.addObjective(n_obj, active=True)


